
source("1.2.1.2.1_Hessian_RD.R")

## Sandwich estimator for variance of RD

IF.alpha.beta = function(y, x, va, vb, vc, alpha.ml, beta.ml) {
    ########################################
    
    ### 1. - E[dS/d(alpha.ml,beta.ml)] ############################## Computing
    ### the Hessian:
    
    Hrd = Hessian2RD(y, x, va, vb, alpha.ml, beta.ml)
    hessian = Hrd$hessian
    p0  = Hrd$p0; p1 = Hrd$p1; pA = Hrd$pA
    s0  = Hrd$s0; s1 = Hrd$s1; sA = Hrd$sA
    rho = Hrd$rho
    dl.by.dpA      = Hrd$dl.by.dpA
    dp0.by.dphi    = Hrd$dp0.by.dphi
    dp0.by.drho    = Hrd$dp0.by.drho
    drho.by.dalpha = Hrd$drho.by.dalpha
    dphi.by.dbeta  = Hrd$dphi.by.dbeta
    dpA.by.drho    = Hrd$dpA.by.drho
    dpA.by.dalpha  = Hrd$dpA.by.dalpha
    dpA.by.dphi    = Hrd$dpA.by.dphi
    dpA.by.dbeta   = Hrd$dpA.by.dbeta
    
    S = cbind(dpA.by.dalpha, dpA.by.dbeta) * dl.by.dpA
    
    IF_alpha_beta = solve(hessian) %*% t(S)
    
    return(IF_alpha_beta)
    
} 
